Roughness induced boundary slip in microchannel flows 
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Surface roughness becomes relevant if typical length scales of the system are comparable to the 
scale of the variations as it is the case in microfluidic setups. Here, an apparent boundary slip is often 
detected which can have its origin in the assumption of perfectly smooth boundaries. We investigate 
the problem by means of lattice Boltzmann (LB) simulations and introduce an "effective no-slip 
plane" at an intermediate position between peaks and valleys of the surface. Our simulations show 
good agreement with analytical results for sinusoidal boundaries, but can be extended to arbitrary 
geometries and experimentally obtained surface data. We find that the detected apparent slip is 
independent of the detailed boundary shape, but only given by the distribution of surface heights. 
Further, we show that the slip diverges as the amplitude of the roughness increases. 

PACS numbers: 83.50.Rp,68.08.-p 
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In microfluidic systems the surface to volume ratio is 
large causing boundary effects to be significantly more 
important than in macroscopic devices. Since even on 
atomic or molecular scales a perfectly smooth surface is 
an idealized model, the shape of the boundary is an im- 
portant property. Additionally, it is of technological in- 
terest to design surfaces with well defined structures and 
properties [ITS] • A commonly investigated surface prop- 
erty is the apparent slip originating for example from the 
surface wettability, electrostatic interactions, impurities, 
or surface structuring [3j. Navier characterized hydro- 
dynamic slip by postulating that the fluid velocity v(x) 
at the boundary (x — 0) is proportional to the shear 
rate ^ and the slip length (3 For macroscopic sys- 
tems the simple no-slip boundary condition (/? = 0) is a 
valid assumption. However, if the height of surface vari- 
ations is not small compared to typical length scales of 
the system, the position of the boundary is not clearly 
defined and experiments might detect slip due to not 
accurately determined wall positions. The influence of 
roughness on the slip length [3 has been investigated by 
numerous authors. Roughness leads to higher drag forces 
and thus to no slip on macroscopic scales, as shown by 
Richardson [f| and Jansons Q . This was experimentally 
demonstrated by McHale and Newton Q. Jabbarzadeh 
et al. performed molecular dynamics (MD) simulations 
of Couette flow between sinusoidal walls and found that 
slip appears for roughness amplitudes smaller than the 
molecular length scale [§]. Also, roughness can cause 
pockets to be filled with vapor or gas nano bubbles lead- 
ing to apparent slip 0, @|- Recently, Sbragaglia et al. 
applied the lattice Boltzmann (LB) method to simulate 
fluids in the vicinity of microstructured hydrophobic sur- 
faces [lOj] and Varnik et al. have shown that even in 
small geometries rough channel surfaces can cause flow 
to become turbulent. A common setup to measure slip 
is to utilize a modified atomic force microscope (AFM) 
to oscillate a colloidal sphere in the vicinity of a bound- 



ary 3 Hi Id | • Vinogradova and Yakubov demonstrated 
that assuming a wrong position of the surface during 
measurements can lead to substantial errors in the deter- 



mined slip lengths LJ| . They showed that measurements 
can be interpreted by assuming a modified sphere radius 
instead of Navier's slip condition, so that the position of 
a no slip wall would be between peaks and valleys of the 
rough surface. In this paper we follow this idea. We an- 
swer the question at which distinct position the "effective 
boundary" has to be placed and study the influence of a 
wrongly determined wall position numerically. 

Panzer et al. gave an analytical equation for f3 for small 
cosine-shaped surface variations [15}. It is applicable to 
two infinite planes separated by a distance 2d being much 
larger than the highest peaks h max . Surface variations 
are determined by peaks of height h max , valleys at h m i n 
and given by h(z) = h max /2 + /i max /2 cos(gz). Here, q 
is the wave number and the corresponding slip length is 
found to be 
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Higher order terms cannot easily be calculated analyt- 
ically and are neglected. Thus, Eq. Q] is valid only for 
k = qh max /2 <C 1. However, for realistic surfaces, k can 
become substantially larger than 1 causing the theoret- 
ical approach to fail. Here, only numerical simulations 
can be applied to describe arbitrary boundaries. 

We use a 3D LB model as presented in 16|, [17, H3 to 
simulate pressure driven flow between two infinite rough 
walls. Previously, we applied the method to study flows 
of simple fluids and complex mixtures containing surfac- 
tant in hydrophobic microchannels [lil Il9j |. Here, we 
only shortly describe our method and refer to the litera- 
ture for details. The lattice Boltzmann equation, ?7i(x + 
Ci,t+l)—r)i(x, t) = fli, with i = 0, 1, . . . , b, describes the 
time evolution of the single-particle distribution rjifat), 
indicating the amount of quasi particles with velocity Cj, 
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at site x on a 3D lattice of coordination number b = 19, at 
time-step t. We choose the Bhatnagar-Gross-Krook col- 
lision operator = — T _1 (?7 i (x, t) — ?7- cq (u(x, t), ry(x, t))), 
with mean collision time r and equilibrium distribution 
Vi q 3 3 • Simulation lattices are 256 lattice units long 



in flow direction and the planes are separated by 62 sites. 
Periodic boundary conditions are imposed in the remain- 
ing direction allowing us to keep the resolution as low as 
16 lattice units. A pressure gradient is obtained as de- 
scribed in [18]. An effective boundary position can be 



found by fitting the parabolic flow profile 
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via the distance 2d = 2d c g. j3 is set to here and vis- 
cosity \i as well as pressure gradient ^ are given by the 
simulation. To obtain an average value for d e ff, a suffi- 
cient number of individual profiles at different positions 
z are taken into account. Alternatively, the mass flow 
/ v (x)pdx can be computed to obtain 2d ff. Both meth- 
ods are equivalent and produce identical results. The 
so found d e g gives the position of the effective bound- 
ary and the effective height h c g of the rough surface is 
then defined by d — d c g (see Fig. [IJ. As rough model 



FIG. 1: (Color online) The effective boundary height /i c ff 
is found between the deepest valley at ft m in and the highest 
peak at /i max . 

surfaces we choose a randomly generated roughness and 
three periodic ones for which the average height or av- 
erage roughness R a is given by h max /2. Cosine-shaped 
boundaries are given by h(x) = h max /2 + h max /2 cos(qx), 
squares have a height of /i max and are separated by /i max 
lattice sites. Triangular structures are 2/i max wide and 
have a height of /i max (see Fig. [2]). Randomly generated 
surface structures are created by choosing for every lat- 
tice position of the boundary the height h{x) as a random 
integer number between and /i max . For determining 
h e s we average over five surfaces generated with differ- 
ent sequences of uniformly distributed random numbers. 
All wall types are geometrically similar, i.e., the effective 
height h e g scales linearly with /i max . 

In Fig. [3J the effective height h c s obtained from our 
simulations is plotted versus R a for cosine shaped sur- 
faces with qh max /2 = k = lj jj> § (symbols). Lines are 
given by the analytical solution of Eq. Q] For k < 1 the 
simulated data agrees within 2.5% with Panzer's predic- 
tion. However, for k — 1 a, substantial deviation between 
numerical and analytical solutions can be observed be- 
cause Eq. [1] is valid for small k only. The inset of Fig. [3] 
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FIG. 2: Periodic surfaces: a) cosines given by h(x) = 
hmax/2 + ftmax/2 cos(qx) . b) squares with height and sepa- 
ration given by /i max . c) triangles, /i max high and 2/i max wide. 



depicts the ratio of /3//i max according to the theory of 
Panzer. In the case of large k > 1, the theory is not able 
to correctly reproduce the increase of (3 with increasing 
fomax anymore. Instead, f3/h ma , x becomes smaller again 
due to missing higher order contributions in Eq. [TJ Our 
simulations do not suffer from such limitations allowing 
us to study arbitrarily complex surface geometries. 
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FIG. 3: (Color online) Effective height h c g over average 
roughness R a for a cosine geometry and different variations k. 
Symbols denote numerical data and lines are given by Eq. [l] 
The inset shows /3(fc)//i max according to equation (fTJ). For 
k > 1 the slope becomes negative, demonstrating that the 
theory fails for more complex surface structures. 

In Fig. h e s is plotted versus R a for different types of 
roughness. By performing a linear fit to the data as given 
by the lines we find for the uniformly distributed rough- 
ness that the position of the effective wall is at c = 1.84 
times the average height of the roughness R a = h max /2 
or at 92% of the maximum height h max . For squares and 
triangular structures we find constants of proportional- 
ity of c = 1.90 and c = 1.69 indicating that the shape of 
the surface variations indeed affects the position of the 
effective boundary. However, the effect of the shape is 
small compared the effect of the height of the variations. 
All surface structures are geometrically similar causing 
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the linear dependence between h e s and R a — h max /2. 
When converting our 3D random roughness into a purely 
2D structure, the difference in the measured constant of 
proportionality c is in the range of the error of the fit 
algorithm. This is a surprising result since in three di- 
mensions the flow can pass sidewise a roughness element. 
The measured h c g is found to be independent of the flow 
velocity over more than three decades and does not de- 
pend on the pressure either, i.e., h c g is independent of 
the Reynolds number. 
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FIG. 4: (Color online) a) Effective height h c s versus R a for 
triangles, blocks (see Fig. [5}, and an equally distributed ran- 
dom roughness, b) h e g versus R a for triangles with /i max = 5 
and 10. The distance between triangles a is varied to obtain 
the given R a . Values of /i ma x = 5 are scaled by a factor of 2. 

In reality high pikes on a smooth surface may occur, 
so that the average roughness R a is much smaller than 
the maximum height /i max . To observe such cases we 
simulate a triangle geometry with additional void space 
a between the roughness elements. As maximum height 
/i max we choose 5 and 10 lattice sites. In similarity to 
Fig. 2] we plot the effective surface height h c e over the 
average roughness R a in Fig. [5J In this case the average 
roughness is smaller than the half of /i max , i.e., R a — 
< The values of /i max = 5 are scaled by 
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a factor of 2 to fit them with the values of /i max = 10. Due 
to the geometrical similarity of the surface structure this 
scaling is possible. For comparison with Fig. |4ji the linear 
fit with slope c = 1.69 is plotted. In Fig.[4jD we see that 
the maximum height has the strongest influence on the 
effective height h g and not the distance a. For small R a 
created by a large additional distance a, h c g converges to 
zero corresponding to a flat surface. For small a the data 
converges to the triangle geometry as given in Fig. 0£t. 
For a medium a sa 2/i max the effective wall is still in the 
range of 75% of the maximum height /i max . This is an 
important result, since it demonstrates that the distance 
between h e g and h m { n can be much larger than R a and 
that in such cases h e g > 6 • R a can be obtained. On the 
other hand, for large a this results in h e g < 0.7 • h max . 
Therefore, in the case of large a rj 2/i max , the effective 
wall h c s cannot be approximated by the maximum height 
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FIG. 5: (Color online) Simulated h c s versus R a for a gold 
coated glass surface and a randomly generated surface with 
Gaussian distributed heights. The background image visual- 
izes the gold surface (left) and the artifically generated struc- 
ture (right). 

We obtained AFM data of a gold coated glass surface 
with a maximum peak to valley distance of 64nm. The 
sample size is l^m 2 represented by 512 x 512 data points. 
A lattice constant of the LB simulation can be scaled to 
1.9nm by setting the relaxation time r to 1.15 and by 
mapping the speed of sound and the viscosity to the val- 
ues for water (c s = 1.5 • 10 3 m/s, fi = 1.02 • 10~ 6 m 2 /s). 
ft e ff can then be measured as in previous paragraphs of 
this paper by loading the AFM data onto our simula- 
tion lattice. For the simulations presented in this para- 
graph, the channel width is set to 128 lattice units. The 
simulated effective height of the gold surface is depicted 
by the square at i? a =21nm in Fig. O Data points at 
R a = 4 and 8 are obtained by downscaling the original 
data set. We find that the distribution of surface heights 
follows a Gaussian distribution and use this distribution 
to generate an artificial random surface with identical 
height distribution. In contrast to the AFM data, our 
data points are fully uncorrelated, while the gold surface 
shows distinct structural properties as can be observed in 
the background images of Fig. [5] For artificial surfaces, 
the average roughness R a can be scaled by scaling the 
width of the distribution of random numbers allowing us 
to determine h e s for R a up to 40nm. As shown by the 
dotted line, the measured h e s linearly depends on R a 
with a constant of proportionality of c = 1.43. The data 
obtained from the gold coated surface follows the same 
linear dependence demonstrating that the actual shape 
of a surface does not influence the effective surface po- 
sition, but only the distribution of heights needs to be 
known. 

The most important question to be answered by our 
simulations is the effect of a wrongly assumed position 
of a surface on experimental measurements. As men- 
tioned in the introduction many groups use an approach- 
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ing method to measure the slip length (3. Here, a colloidal 
sphere at the tip of a cantilever immersed in a fluid is os- 
cillated in the vicinity of a surface, or the two cylinders 
of a surface force apparatus (SFA) are brought close to 
each other. The distance between the surfaces can be- 
come very small - even down to contact. To study the 
influence of the roughness on an apparent slip effect, we 
assume the surface to be placed at h max as it is com- 
monly done in experiments (l3| . Then, we measure the 
slip length (3 by fitting Eq. [2] The wrong position of 
the surface causes a substantial error in the detected slip 
as can be inferred from Fig. [6l Here, j3 is given versus 
R a for randomly generated boundaries with the heights 
of the surface obstacles following the Gaussian distribu- 
tion given by the AFM data of the gold surface. For 
small R a (and thus large separation of the plates) (3 is 
in the range of h max — h e s and can be neglected in most 
practical cases. However, the detected slip diverges if R a 
becomes large and grows to 80nm for R a = 55nm. Here, 
a large R a is equivalent to the channel width becoming 
very small - an effect also common in typical surface ap- 
proaching experiments or microchannel flows. For curved 
surfaces, as they are utilized in surface force apparatuses 
or AFM based slip measurements, the detected (3 can be 
even larger due to higher order components of the flow 
field. This might explain experiments reporting large slip 
lengths of (3 rj lOOnm 0, H|. 
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FIG. 6: (Color online) Slip length /3 versus R a for water and 
the randomly distributed roughness. The line is a guide to the 
eye. By assuming h a g — /i ma x as it is common in experiments, 
f3 is in the range of /i ma x — h g for small R a , but for larger 
R a the apparent slip diverges. 

In conclusion we performed LB simulations of pres- 
sure driven flow between two rough plates. By varying 
the roughness we found that there exists an imaginary 
effective plane where the no slip boundary condition is 
valid. We compared our results to analytic calculations 
of Panzer et al. and found good agreement in the case of 
small variations (k < I). Large and more realistic per- 
turbations (k > I) can only be covered by simulations 
as presented in this paper. By simulating flow of water 
along a gold coated surface and a randomly generated 
one with identical height distribution, we demonstrated 
that the position of the effective plane is independent of 
the actual boundary structure and that only the distribu- 
tion of heights is relevant. We showed that apparent slip 



due to errornous assumptions of the surface structure can 
become very large if the distance between the boundaries 
is small - as it is typical in dynamic microfluidic exper- 
iments. Our simulations can be of practical importance 
for experimental measurements of boundary slip induced 
for example by electrostatic interactions, surface wetta- 
bility or impurities. Due to the precise measurements 
needed, ignoring the influence of surface roughness leads 
to substantial errors in the determined slip. A simulation 
of the flow along a surface generated from AFM data al- 
lows to determine how an experimentally detected slip 
might have to be corrected in order to take the surface 
structure into account. 
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